function gridStructure = constructTriangleGrid(sites,edges,triangles)

Nx = 30; Ny = 30;
minx = min(sites(:,1)); maxx = max(sites(:,1));
miny = min(sites(:,2)); maxy = max(sites(:,2));

dx = (maxx-minx+1)/Nx;
dy = (maxy-miny+1)/Ny;
gridStructure = repmat(struct('triangles',[], 'edges',[]),Nx,Ny);

% insert edges
for i=1:size(edges,1)

    ax = sites(edges(i,1)+1,1); ay = sites(edges(i,1)+1,2);
    bx = sites(edges(i,2)+1,1); by = sites(edges(i,2)+1,2);
    
    % add edge in grid where the ends lie
    indax = floor((ax-minx)/dx) + 1; inday = floor((ay-miny)/dy) + 1;
    indbx = floor((bx-minx)/dx) + 1; indby = floor((by-miny)/dy) + 1;
    
    if(isempty(find(gridStructure(indax,inday).edges==i)))
        gridStructure(indax,inday).edges = [gridStructure(indax,inday).edges; i]; 
    end    
    if(isempty(find(gridStructure(indbx,indby).edges==i)))
        gridStructure(indbx,indby).edges = [gridStructure(indbx,indby).edges; i];
    end
    
    % check whether edge crosses another grid
    slope = (by-ay)/(bx-ax);
    for j=min(indax,indbx)+1:max(indax,indbx)
        indy = slope*(((j-1)*dx+minx)-ax) + ay;
        k = floor((indy-miny)/dy) + 1;
        if(isempty(find(gridStructure(j,k).edges==i)))
            gridStructure(j,k).edges = [gridStructure(j,k).edges; i];
        end
    end
end

% insert triangles
for i=1:Nx
    i
    for j=1:Ny          
        for k=1:length(gridStructure(i,j).edges)
            searchEdge = gridStructure(i,j).edges(k);
            searchVertex1 = edges(searchEdge,1)+1;
            searchVertex2 = edges(searchEdge,2)+1;
            
            for l=1:length(triangles)
               if(~isempty(find(triangles(l,:)==searchVertex1)) || ~isempty(find(triangles(l,:)==searchVertex2)))
                   if(isempty(find(gridStructure(i,j).triangles==l)))
                       gridStructure(i,j).triangles = [gridStructure(i,j).triangles; l];
                   end
               end
            end
        end
    end
end

end